{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 10.6变分Logistic 回归"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "选择$p(w)=\\mathcal N(w|m_0,S_0)$的高斯共轭先验，假定超参数已知"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "边缘似然函数:$$p(t)=\\int \\bigg [\\prod_{n=1}^Np(t_n|w)\\bigg ]p(w)dw$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "利用Logistic sigmoid函数的变分下界，t的条件概率可以写成\n",
    "$$p(t|w)=e^{at}\\sigma (-a)\\ge e^{at}a(\\xi)exp\\bigg\\{-\\frac{a+\\xi}{2}-\\lambda(\\xi)(a^2-\\xi^2)\\bigg \\}$$\n",
    "\n",
    "\n",
    "因此存在⼀个变分参数ξn，对应于训练集的每个观测$(\\phi_n; t_n)$。使⽤$a =w^T\\phi $，乘以先验概率分布，我们可以得到下⾯的t和w的联合概\n",
    "率分布。\n",
    "\n",
    "$$p(t,w)=p(t|w)p(w)\\ge h(w,\\xi)p(w)$$\n",
    "\n",
    "\n",
    "$$h(w,\\xi)=\\prod_{n=1}^N\\sigma(\\xi_n)exp\\{w^T\\phi_nt_n-(w^T\\phi_n+\\xi_n)/2-\\lambda(\\xi_n)([w^T\\phi_n]^2-\\xi^2_b)\\}$$\n",
    "\n",
    "代入$p(w)$得到变分后验概率分布的近似\n",
    "$$q(w)=\\mathcal(w|m_N,S_N)$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10.6.2最优化变分参数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$lnp(t)=ln\\int p(t|w)p(w)dw \\ge ln \\int h(w,\\xi)p(w)dw=\\mathcal L(\\xi)$$\n",
    "\n",
    "将$w$看作一个潜在变量，使用EM算法，用$\\xi^{old}$计算w上的后验概率分布，在最大化完整数据似然函数的期望\n",
    "$$Q(\\xi,\\xi^{old}）=\\mathbb E[ln\\{h(w,\\xi)p(w)\\}]$$\n",
    "代入$h(w,\\xi)$:\n",
    "$$Q(\\xi,\\xi^{old}=\\sum_{n=1}^N\\bigg \\{ln\\sigma(\\xi_n)-\\frac{\\xi_n}{2}-\\lambda(\\xi)(\\phi^T\\mathbb E[ww^T]\\phi_n-\\xi^2_n\\bigg\\}+const$$\n",
    "\n",
    "对$\\xi_n$求导再代入$\\sigma(\\xi)$和$\\lambda(\\xi)=\\frac{1}{2\\xi}[\\sigma(\\xi)-\\frac{1}{2}]$可求出$\\xi^{new}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "计算如下:\n",
    "* 初始化变分参数$\\xi^{old}$\n",
    "* E-step:$$m_N=S_N(S_0^{-1}m_0+\\sum_{n=1}^N(t_n-\\frac{1}{2})\\phi_n)$$\n",
    "$$S_N^{-1}=S_0^{-1}+2\\sum_{n=1}^N\\lambda(\\xi_n)\\phi_n\\phi_n^T$$\n",
    "* M-step:\n",
    "$$(\\xi_n^{new})^2=\\phi_n^T\\mathbb E[ww^T]\\phi_n=\\phi_n^T(S_N+m_Nm_N^T)\\phi_n$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10.6.3超参数的推断"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$p(w|\\alpha)-\\mathcal N(w|0,\\alpha^{-1}I)$$\n",
    "$$p(\\alpha)=Gam(\\alpha|a_0,b_0)$$\n",
    "边缘似然函数变为:\n",
    "$$\\begin{align}p(t)&=\\iint p(w,\\alpha,t)dwd\\alpha\\\\&=\\iint p(t|w)p(w|\\alpha)p(\\alpha)dwd\\alpha\\end{align}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "引入变分分布$q(w,\\alpha)$:\n",
    "\n",
    "$$\\mathcal L(q)=\\iint q(w,\\alpha)ln \\frac{p(w,\\alpha,t)}{q(w,\\alpha)}dwd\\alpha$$\n",
    "$$KL(q\\lVert p)=-\\iint q(w,\\alpha)ln(\\frac{p(w,\\alpha|t)}{q(w,\\alpha)}dwd\\alpha$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "继续利用Logistic sigmoid函数的变分下界:\n",
    "$$lnp(t)\\ge\\mathcal L(q)\\ge\\tilde L(q,\\xi)=\\iint q(w,\\alpha)ln\\bigg\\{\\frac{h(w,\\xi)p(w,\\alpha)p(\\alpha)}{q(w,\\alpha)}\\bigg\\}dwd\\alpha $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "假设变分分布可在参数和超参数上分解\n",
    "$$\\begin{align}lnq(w)&=\\mathbb E_{\\alpha}[ln\\{h(w,\\xi)p(w|\\alpha)p(\\alpha)\\}]+const\\\\&=lnh(w,\\xi)+\\mathbb E_{\\alpha}[lnp(w|\\alpha)]+const\\end{align}$$\n",
    "然后代入$h(w,\\xi)$,$p(w|\\alpha)$后得到关于w的二次函数，可看作高斯分布\n",
    "$$q(w)=\\mathcal N(w|\\mu_N,\\sigma_N)$$\n",
    "再次代入$p(\\alpha)$\n",
    "$$lnq(\\alpha)=\\mathbb E_{w}[lnp(w|\\alpha)]+lnp(\\alpha)+const$$\n",
    "结果为一个Gamma分布"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "像上节一样使用EM算法更新参数\n",
    "*   初始化变分参数$\\xi^{old}$\n",
    "*   E-step:$$\\Sigma^{-1}_N\\mu_N=\\sum_{n=1}^N(t_n-\\frac{1}{2})\\phi_n$$\n",
    "$$\\Sigma^{-1}_N=\\mathbb E[\\alpha]I+2\\sum_{n=1}^N\\lambda(\\xi^{old})\\phi_n\\phi_n^T$$\n",
    "$$a_N=a_0+\\frac{M}{2}$$\n",
    "$$b_N=b_0+\\frac{1}{2}\\mathbb E_w[w^Tw]$$\n",
    "    其中: $$\\mathbb E[\\alpha]=\\frac{a_N}{b_N}$$\n",
    "$$\\mathbb E_w[w^Tw]=\\Sigma_N+\\mu\\mu_N^T$$ \n",
    "         \n",
    "*    M-step:$$(\\xi_n^{new})^2=\\phi^T_n(\\Sigma_N+\\mu_N\\mu_N^T)\\phi_n$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "#sys.path\n",
    "sys.path.append(\"../\")\n",
    "import matplotlib.animation as animation\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "\n",
    "from prml.rv import VariationalGaussianMixture\n",
    "from prml.features import PolynomialFeatures\n",
    "from prml.linear import (\n",
    "    VariationalLinearRegressor,\n",
    "    VariationalLogisticRegressor\n",
    ")\n",
    "\n",
    "np.random.seed(666)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_toy_data(add_outliers=False, add_class=False):\n",
    "    x0 = np.random.normal(size=50).reshape(-1, 2) - 3.\n",
    "    x1 = np.random.normal(size=50).reshape(-1, 2) + 3.\n",
    "    return np.concatenate([x0, x1]), np.concatenate([np.zeros(25), np.ones(25)]).astype(np.int)\n",
    "x_train, y_train = create_toy_data()\n",
    "x0, x1 = np.meshgrid(np.linspace(-7, 7, 100), np.linspace(-7, 7, 100))\n",
    "x = np.array([x0, x1]).reshape(2, -1).T\n",
    "feature = PolynomialFeatures(degree=1)\n",
    "X_train = feature.transform(x_train)\n",
    "X = feature.transform(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAFeCAYAAACl2PUiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XecVPW9//HXZ6fsLixVBJEiCEgVERFr7BqMLT2aa4ma\nGE1Mcn9Jbq7G9HZNbprGJEqiaZbEa4nGYI0x2AEV6SggCkSKSGd3p31+f5yZZXbK7uzO7Jwyn+fj\nsQ93p5zzXYT3fvZzvt/vEVXFGGOMf9S5PQBjjDFdY8FtjDE+Y8FtjDE+Y8FtjDE+Y8FtjDE+Y8Ft\njDE+Y8FtjDFlEpHbRGSziCwp8ryIyI0iskpEFonI9KznZonIyvRz15RyPgtuY4wp3++BWR08fyYw\nLv1xBfBrABEJAb9MPz8JuEBEJnV2MgtuY4wpk6rOBd7t4CXnAX9UxwtAfxEZCswEVqnqGlWNAX9O\nv7ZDFtzGGNPzhgHrsr5en36s2OMdCld0aMYY40HvPeNo3bp1R7ff/9LLK5YCLVkPzVbV2WUPrJss\nuI0xgbd16w7mPfe7br8/1HBMi6rOKGMIG4ARWV8PTz8WKfJ4h6xVYowxPe9B4OL07JKjgR2q+jYw\nHxgnIqNFJAqcn35th6ziNsaYMonIXcBJwCARWQ98E6eaRlVvBuYA7wNWAXuBS9PPJUTkauBRIATc\npqpLOzufBbcxxpRJVS/o5HkFPlvkuTk4wV4ya5UYY4zPWHAbY4zPWHAbY4zPWHAbY4zPWHAbY4zP\nWHAbY4zPWHAbY4zPWHAbY4zPWHAbY4zPWHAbY4zPWHAbY4zPWHAbY4zPWHAbY4zPWHAbY4zPWHAb\nY4zPWHAbY4zPWHAbY4zPWHAbY4zPWHAbY4zPWHAbY4zPWHAbY4zPWHAbY4zPWHAbY4zPWHAbY4zP\nhN04af/9BurQEcPdODVJbQGgIRRx5fzdonsBCNfVuzwQl2gLddLo9iiMS156ecU7qrq/2+PwEleC\ne+iI4fzh8QfdODUA22MrAZjQ/wDXxtBVqfgiBjeMcXsYrpD4cuqjk90ehnFJqOGYN90eg9fUZKuk\nf3Q8ACu2b3R5JF2zuWW120NwhUYm0hpb6vYwjPGMmgxu2BfeflEXmer2EIwxHlGzwQ1OeFvVbYzx\nm5oO7gy/hLdV3cYYsOD2Xb+7LjK1Jqtu63Mbs0/NBzf4L7whOC0T1WZU97g9DGN8xYI7zU/hHYSW\niaa2ojt/ANuuhG2fQXd8HU1ucHtYxnSLiMwSkZUiskpErinw/AARuV9EFonIPBGZkvXcF0RkiYgs\nFZH/LOV8FtxZ/DbTxK9Vt5KEnd+BxEog6Xwk18LO76DpxUbG+IWIhIBfAmcCk4ALRGRSzsu+CixU\n1anAxcAN6fdOAT4FzAQOA84WkbGdndOCO4dfZpr4uuqOvZpeDZrMelBBE9D6XNG3WZ/beNRMYJWq\nrlHVGPBn4Lyc10wCngRQ1RXAKBEZAkwEXlTVvaqaAP4FfLCzE7qyctIPVmzf6IuVlZtbVvtvRWVq\nixPSeVoh5f0fmsZ/Utpc7g/9QSKyIOvr2ao6O/35MGBd1nPrgaNy3v8qTiA/LSIzgYOA4cAS4Psi\nsh/QDLwPWEAnLLgL6B8dz/bYSs+Hd11kKqn4IreH0XXhg0BCoPH2j0sDhH32Q8j4gzSgkYnlHOEd\nVZ1RxvuvB24QkYXAYuAVIKmqy0Xkh8BjwB5gIe1/FS3IWiVF+Olipe963eHxEBoJRLMfhLr+EO38\n34a1S4zHbABGZH09PP1YG1XdqaqXquo0nB73/sCa9HO3quoRqnoCsA14rbMTWnB3wA8XK/3Y6xYE\n+l4DjbOgbgBIX6g/Gfp+C6HjXRvLrJqM6QnzgXEiMlpEosD5QLtd9ESkf/o5gE8Cc1V1Z/q5wen/\njsRpp9zZ2Qkr0ioRkf7Ab4EpgAKXqerzlTi225yLlSs93zLZ3OKv3QOFKDR+1PkwxsdUNSEiVwOP\nAiHgNlVdKiJXpp+/Geci5B9ERIGlwOVZh7g33eOOA59V1e2dnbNSPe4bgEdU9cPpnyq9KnRcz/B6\nvxt8eqHSmABQ1TnAnJzHbs76/HngkCLvfU9Xz1d2q0RE+gEnALemBxEr5SeGn/ih3+3Hlkk5rM9t\nalkletyjgS3A70TkFRH5rYj0rsBxPcUP4Q0+vFDZDdbnNrWuEsEdBqYDv1bVw3GmtBRa8nmFiCwQ\nkQXbt26twGmrz+vhXWtVtzG1qhLBvR5Yr6ovpr++ByfI21HV2ao6Q1Vn9N9vvwqc1h1+mGlSC1U3\nWLvE1K6yg1tVNwLrRCSTaKcCy8o9rpd5eVl8rVTd1i4xtaxS87g/B9whIouAacAPKnRcT/NyeNdK\n1W1MLapIcKvqwnQbZKqqvl9Vt1XiuF7m9X431E7LxJhaYysny+Dl8K6FlontFmhqlQV3mbwc3mBV\ntzFBZMFdAV6daVILVbcxtciCu0K8PNMkyFW3tUtMLXIluFti8c5f5FNeC2+ruo0JHleCuyEUYcWG\nzazYsNmN0/cYr/a7bXqgMcHiWqtkYtNQAAvvKgpqeFu7xNQaV3vcFt7VYy0TY4LD9YuTQQ9vrwlq\n1Q22d4mpHa4HNzjhPbFpaOD63l6baRLkqtv2LjG1xBPBnRHU6ttL4Q3BrrqNqQWeCm4IXnh7rd+d\nqbqDGt7WLjG1wHPBDRbePS2oLRNrl5ha4cnghuD1vb0W3hDcqtuYoPNscGcEqfr2UngHueq2dokJ\nOs8HNwQzvL3Cqu4yJNdAy6MQmwck3BuH7gLd7d75TdWF3R5AqSY2DWX57rfbwnvCsMEuj6j7nGmC\nK5nQ/wBXx1EXmUoqvojNLasZ3DDG1bFUWmtsKfXRyT109CS6++cQXw4kgTDsiSB9vwahA3vonIWG\nsQ7d/WtIve18HRqN9L4KQvtXbwzGFb6ouDMyfW8IRvVtLZOe0dMXKbXlCYgvA1pxKu0WYDe6+xc9\net72g9iD7vwepNalx5CA5Cp013dxfpiYIPNVcGcEIby91O+GYLZMyup1azMk1zptiLwDPwXEct8A\nqU2Qeqf75+wCbX2O/IBWZ9yxhVUZg3GPL4MbLLwryarunPc2/x+6/Wp05w/Q7Z9Hd/8KyN6KuFhF\nK6BZz+lupw/eE/3n1Gacij9Xomo/PIx7fBvc0D68/RrgXgrvIFbdXaUt/3QuOBIDmoEExBege27f\n96LocUAk/83SF0JDgBS653dO6O+83vnvnlupZAtDwmOB+gLPhCA8umLnMd7k6+CGYPS9vTTTJEjh\n3a2pga0PkV/JxiH2NJmZI9J4JtQNZ19wRoEGpOlq57zND0DsGed9NKff/xzafH83v5MCokdA3X60\nn18QhdBoCB9SufMYT/J9cGcEIby9UHXXvFSxtoaCtqQ/jyL9vunM4Kg/Exo/gvT/KYTTM3NaHiO/\nBx6DlscrONAw0vdbUH86SH+Q/aDhbKTvVyp4DuNVgQlu8H94g/stE6jxqjs8DpD8x6UfSFPWAyGI\nHoH0+jjSMAukT9Zze4scvLn0cZRCGp3z9/8F0v/nSOMHKNjCMYETqOAGfy+V90K/O+ibUHVGel2A\n0/rIhLcAUaTXJ0o/SGhU4cfrRpYzNGPaBC64M/xafXspvIOiS1V3aBjS9/sQPR7qDoTINKTPVyE6\nreTzSa+Lcfrf2eFfj/S+uIsjLyK1Iz1VsdCsEuMGEZklIitFZJWIXFPg+QEicr+ILBKReSIypdT3\nFuKblZPdkb3a0k8rLftHx7M9tpIV2ze6uroyiCsqSxIagvS+ovvvD49B+n4bbX4Qkm9CaCTSeC6E\nhpc3Lm1F9/wa4q/itESS0HCec2zjGhEJAb8ETgfWA/NF5EFVXZb1sq8CC1X1AyIyIf36U0t8b57A\nVtwZVnl3T01X3ZUQGoY0XYX0ux5p+kz5oQ3ont+kQzuB0y+PQcuDEHuh7GObsswEVqnqGlWNAX8G\nzst5zSTgSQBVXQGMEpEhJb43T6Ar7oxMeC/f4Ozp4JfqO1N5u8WZ270oUFV3z+5h0oO0GeIvkb+Z\nVSva/BASPdqNUflGItVa7nWbQSKyIOvr2ao6O/35MGBd1nPrgaNy3v8q8EHgaRGZCRwEDC/xvXlq\nIrgz/Ng68cKGVEFpmWhkIhJf7vYwukf3UPQXZN1R1aH4kvQq97fId1R1Rhnvvx64QUQWAouBVyhj\nRVbFWiUiEhKRV0TkoUodsyf4sXXi5hzvoLVMwGO3N9MWaH0K3Xs7tD5N/vzvtLqBFF4pKRCe1IMD\nNCXYAIzI+np4+rE2qrpTVS9V1WnAxcD+wJpS3ltIJXvcXwB8Uc74MbzB3X53UKYHeur2Zskt6I4v\npUP7UXTvH9DtX4bUtgIvroPeF+FMVcwIAY1I44eqM15TzHxgnIiMFpEocD7wYPYLRKR/+jmATwJz\nVXVnKe8tpCLBLSLDgbOA31bieNXgt/nebl+shGDN7fZC1a17f5/efTAzra8VdIcT5AVI9Bikz1cg\nfJgzVTF6EtLv+xDyR9svqFQ1AVwNPIpTvN6tqktF5EoRuTL9sonAEhFZCZyJU+gWfW9n56xUj/vn\nwFeAPsVeICJXAFcAHHDgsAqdtnx+6nu7OU0wc9OFIMj0ut29UKmQWOz8t50UxF8p/rbweKSPd/a2\nMQ5VnQPMyXns5qzPnwcKbiJT6L2dKbviFpGzgc2q+lJHr1PV2ao6Q1Vn9B84sNzTVpSfWiduVt7W\nMqkkofg/v8DP0jVlqsTfkOOAc0VkLc4cxFNEpPDveh7mpy1i3W6bBCW8weWWSWQGTp86WxiixxR4\ncQKSW2y1pAEqENyqeq2qDlfVUTiN9SdV9cKyR+YCP20R69ZWsEGaZZKpuns8vJNvOVu6Nj/khG+a\n9L4E6gYDDThdy3qoG5reLyVrnM1z0G1XoTuvRbd/Bt37J+z2ZLWtpuZxl8ovfW+35ngHaWFOT8/t\n1r13QesTQBwlBC33QsMHkcZzQPog/X4I8SWQfBtCwyAyiezdCbX1GWi5j3Z7hLc+hRJFen2sx8Zt\nvK2izTRVfUpVz67kMd3ip8rbWibl65GqO7E6HdoxnIuQ6Zv6ttyN7rgOUu8CApFDITojPVc7R8tf\nyb+xQwxaH8eq7tplFXcH/LRUvtozTWpvlkkCbX4MYk8BKYgehzS8D6TQopj0cWMv0v5elVlS69Bd\nP4Cma2H3jem7tdeBNCC9r4RIevO41Pai40FjII0lfY8mWOzydQm8Xn27dbEyiLNMilXeuuvHTpsj\n9bZzN/eWv6E7vwekOjhqR/+81AnlXd+G1Bs4AZ+ex73755BM/10LFbl/pPSx0K5hFtwlsvAuLmjh\nnSfxOiRW0X45etwJ8VjxOddSfwyd3pFGd5M/lzuJtv7DOUav82l/Ywecr3td1PFxTaBZcHeBhXe+\noN0xp+D2r4nVFO4nt6KJVcUPFjoIGs4mf8pfRpzC/wQTkFwHybcgPBrp+w1ntaT0h9B4pOmLSHRm\nKd+OCSjrcXeR1/vebqyuDFK/O6Ndv7tuAM4/ldwtVaNIoQuKWaTxA2jkCNj1Pzj3osy0VuohciTE\nXyz0LkgsQ3d+BySK9P4s0udL5Xw7JmCs4u4mL1ffbrVNglR1Q1a/Ozodp+WRexPhENQf2+nxJDwS\n6f8TaDgH6oZCaAzS+zKk6dMQPY78Xf8Up8JvBd2F7v5ZegaKMQ4L7jL4IbyrJYgtE8iEdwTp+3Wo\nG44T4BGoG4L0vRakd2kHlF5I44eRfj9C+n4Lok7gS+/LnHtU1h3ktEIKtlWSaOvc8r8pExgW3GXy\nenjXcr9bSaK6E+3mfOd24R0aivT7AdLvx0749vtx8RkfXSJQfwLS73tI4wcoHNyJIlu9mlplPe4K\n8HLfu9qrK73Q71YUmh92Fq9oHCSMNpwLjWcjee2OTo6VO8e7k552WcLFfkuqRyI+vN2a6TFWcVeQ\nl6vvmup3t/wDmu8F3QvEnfs1Nv8VWh7r1uE6nuOdomIrGEPDnAuW7XreUagbBtEjKnMOEwgW3BXm\nxfCu9sVK11smLQ+Qv0y81bkjejcoraTCo5yjZMJbd6O7b0K3XYZuuxTd+X1Ilv/nK02fdjafCo1z\n+t4NH0L6XUfxKYWmFrnSKmltyZ1WFSzZm1SBN1on1Z4mmGmZuHKj4VSRm+emdqBoye0STW2HPbMh\n7oR1KjQSaTiNlphS33wHpDbSVm0nV6I7v+XMHil2wTK5Dt17h7OgR3pB/Syk8X20n60iEH0PEn1P\nSWM0tcm1invV65tY9fomt07f47y4RaxblXfVhYYWefyA0kObFOz6bjq0k85Hci26905IvE5r8m3a\nz+tWII62Pl34gMnN6M5vQ2IpEAPdDi33o3v+UPK3ZUyGK8HdGAkzacj+QG0EONRueIMLLZNeH6f9\nTXVxvm78OJr8N9p8F7rnFjQ2r/iMk/hSSO2kff9aQeNo8i1UhdZUktZU9vMxSK4veDhteYj8Dada\nITY3fd/JAlLb0T2znb24t3/e2dM7bxGQqUWu9rgnDdm/JgK8lsO72v1u1RaIvQpSh9OCCEFoJPT5\nIhCDHV+D5keg9WmnDbLz+2ihHfxSm0ELhXoMSIDs22FkX3jXF58imFhD4Q2pIoV749qM7vw6xJ4D\ndoNug5aH0N03dvTtmxrhiYuTtRDgXrs1WhDDW1HYdT20/hO0BSdaFXQHhEbAnt/iBG86aLUFkm9C\n6zP5BwsdRP5KSUAaIDzdOR6RrPBWZ0vWYispQ8MLH48Y1O2f/720PuPMhiGnoo8vLVrVm9rhieDO\nCHqAe63vHbjwTqx2NmdqV0GnINUCzX+jYHBqK7Q+n/94eAyER9O+5RICaYKGY6DPNVB/MkgvlCga\nnUFr43/QGl9TcGjSeBb5OwVGIXIE1PUv8L28Rv7MGJzvIflmwXOY2uGp4M6ohQCH2g3vHpNcl79D\nKuDs+dHB36G6hryHBIG+X4HGWVDXzwns+vdAv+8gRBFpQHpfjAyYjQy8DWn6PESPcs5WaL53aATS\n57+g7kCcHyBRqD8Jabqy8JhCB1J00leBCt3UFk+vnMyE97JNW9rCe+y4IW4OqWK8dF/Lak8V7LEp\ngqEDQKRAeEchNAkSb6bbD9nqIXpywcMJUWj8qPNRotzFOu3uqBOe4NxjkjjOP73iM1yk/mS0ZQ7t\nL0aGnNAOH1LyeEwwebLizhXUCjzTOvFC37talXePtkzCE9LVaHY9IiBRaDgB+nxp351jpAGIQsNp\nEJ1W8aFkB3h+BV5op8Ecdf2Rvl+FuhE4i2/CEDnUeczUPE9X3LmCWoF7pfrOVN49racW5wiC9r0O\n9vweYvMBdarT3pcj0htCvdEBN0JsCbAbwhORukEVO3+uTHhn9joBOrinZQGh0Ui/HzhL9yVM/hRH\nU6t8FdwZQQxwL4V3NTal6rHwliZoutpZQIMiOUvFhQhEDy/5eJp8C/b8Mb3asQEaToXGDyJd+KdT\nNMC1GWIvoKmtSGgsRKdS8Jdg6VXyuUxt8EWrpJigtVC8ctGyWtvB9mTbRKjLC+2u0uQW2PldSKwA\nkqB7nDngu2/p3vEiE9tCvKXln7Rs+6yzBL7lAXTPTeiObzqzXIzphK+DOyNIAe6Vvne1w9uTWh92\ntoVtJwaxl9DU1m4fViMT0ZaHQVtoTe1NL+BphdQGtOVvZQ3Z1IZABHdG0AIc3K2+qxneXrn5QjuJ\ntRRcYi5hSL7d7cNqahskN2WWBwGkl8+30NryZLePa2pHoII7IygBbuHtstBBFNxOVRPO1MNua//P\nTrM+oNhMFGP2CWRwZwQhwGspvME7tz0DoGEWSIHVjtHDujwbRbMml0tdP+emCXlTAqNo/al5Uwkt\nxE2uQAd3ht8D3Ct9b1/P8e4GCQ2Bvl91lr8j6Vklp0DTZ0o+hrbORbd9Ft69CN3+uX03/W36zL45\n5YTSe6CMgsaznfelL2RaiPuDiMwSkZUiskpErinw/H+JyML0xxIRSYrIwPRza0Vkcfq5BSWdT7Xg\nGuEeNXnSVL3zdvcuwizbtKXtc79NI1y+2737WmbmePf0VMHMPSurfgOGDigppIt1ju75A7Q+nvNo\nFHp/Aqk/AU1f6CS5FcIHQ2RiSfuFS3x5u6+7NDfch0INx7ykqjPKOcbU6ZP1oafv7Pb7D2qaVnQM\nIhICXgNOB9YD84ELVHVZkdefA/w/VT0l/fVaYIaqvlPqeMqexy0iI4A/AkNw2nSzVfWGco/bk3Ln\ngfspvN2c712tpfGu3j2niC6HdnJdgdAGiEHzPc6d3YlC9JgujyVThUP7ueEZQQ/y7mhJxnvyN8aZ\nwCpVXQMgIn8GzgMKBjdwAXBXOSesxAKcBPAlVX1ZRPoAL4nI48V+2nhJW4D7bBGPm7dGq+Xw7pKW\nJ4o/l3q3S7dQ60h2iGcUaqfUepiHpKFtW4duGpTTxpitqrPTnw8D1mU9tx44qtBBRKQXMAu4Outh\nBZ4QkSRwS9Zxiyo7uFX1beDt9Oe7RGQ5zjfi+eDO8GOAZy5aulF9W3iXQLcXf076VCS0i546J8wL\nVeVgYd5F75Tbrkk7B3hWVd/Neux4Vd0gIoOBx0VkharO7eggFb04KSKjgMOBFws8d4WILBCRBdu2\nvZv7tCf48QKmW7NO+kfHV2W2idcuWJYsPI2idVHDB6s6lOwLndkf2Rc8cz9Ml2wARmR9PTz9WCHn\nk9MmUdUN6f9uBu7Hab10qGLBLSJNwL3Af6rqztznVXW2qs5Q1RkDBgys1Gkrzo8zUNycMmjhXUTD\ncem53tnzwAXC05HG090aVTvFAr2zULdgzzMfGCcio0UkihPOD+a+SET6AScCD2Q91jvdYkZEegNn\nAEs6O2FFNpkSkQhOaN+hqvdV4phu89tGVm2tkw3Vn3VSjY2p/NY2EaJo329By5MQewHqekP96V3a\n4MpNhXrn2boT3kFtzahqQkSuBh7F+Ul9m6ouFZEr08/fnH7pB4DHVHVP1tuHAPeLCDh5fKeqPtLZ\nOcueDijOGf8AvKuq/1nKe9yeDtgdfppC6NaUwe2xnt9V0ItTBU1pcqcxlqpX30+WPR1w4rSp+ofH\n84rgkh01eHTZY6ikSlTcxwEXAYtFZGH6sa+q6pwKHNsz/FSBuzVlMFN5Q8/N9d7XNllk4e0znVXx\npnRl97hV9RlVFVWdqqrT0h+BCu1suT1wr3LzoiVUZ5Xl5pbV/up7G1MhNbHkvSdkAtzLFzDdWiof\niNugGeNhFtxl8kuAQ3WrbwtvY3qOBXeFeH0KoYV3z9HkJnTPb9AdX0Z3/RBNdO8inDGlsuCuIK/P\nAbfwrjxNboSdX4PWZyC5EeKLYef/orHnqzoOU1ssuHuAlwPcjb53NcO76hctm/8vfZ/IZNaDMdjz\nx/QNi42pPAvuHuTlGSjVrr6rFd5Q5eo7sRIKBbS2Qmpbz5/f1CQL7irw6gVMC+8KkL5FnlCQXj17\nblOzLLiryIvtk+zwrkaAuxXePRbgjecC9TkPRiE6E5HGnjmnqXmuBHdrc4zVS9axesm6zl8cMF7s\nf2f63lCd6jt7Z8Fq9b2hZ6pviR4Nje8H6tO3IYtAdBo0XVbxcxmT4Upw19eHGTvGWYqdCfBaC3Gv\nBjhY66SrpPEcGPgr6Pt1GHAD0vR5JK8KN6ZyXG2VjB0zuO0DajPEvXYBM+jh3VOzToR6JDQSKdrz\nNqZyKrKtayVkwjtjVU54j5kygiCbNGR/z2xgVe0tYqt1R50Mv20Ra0wuz16crMVq3Gvtk2r3vaE6\nlTe4NOfbmArxbHBnq7UQ91L7xK3wrmaAg+11YvzFM62SUmW3VFat3twuvIPWTvHKTYyrub93Jryr\n3ToBZ49vsJs0GO/zRcVdTK1U4l5on1R7qXy1WydQpXnfxlSAr4M7W9BD3CvtkyD3vaHn530bUwm+\na5WUIsjtFC+0T6rdOqlm2yTD2ifGywJTcRcT1Erc7fZJUFda5rL2ifGiQFbcxQStEs+9gXG1q++g\nz/fOyIR3Zu43WAVu3BX4iruYIFXibu8+GPS+d4b1v41X1GxwZysW4n7jZvukVlonYIt3jPssuHNk\nh7gfq3A3Z58EeZ+TQizAjVssuDvg51aKW9V3Lcz3zmUBbqrNgrsEfm2l1Er17XbrJMMC3FSLBXcX\n+bGV4mb1DbXTOsmwADc9zYK7DH5qpbhVfbt5azSvBbiFuKmUmprH3VOKzQ/34txwN1Zets33DvBG\nVR2xeeCm0qzirjC/9MLdrr6rwSutk4zceeBWgZvuqkhwi8gsEVkpIqtE5JpKHNPv/NALd2Phjhvh\n7YULl9kswIOnswwUkZNEZIeILEx/fKPU9xZSdnCLSAj4JXAmMAm4QEQmlXvcIPF6FV7t6rvaUwbB\ne9U37Atw64P7Wxcy8GlVnZb++E4X39tOJSrumcAqVV2jqjHgz8B5FThu4Hg5wN2YeeJm68RLAQ5W\nhftcORnYrfdW4uLkMCA7gdYDR1XguIHV7mKmhza6cmPTqmpuEQveu3CZq9CFTLCLmeVqicV7skAo\nNQOPFZFFwAbgy6q6tAvvbadqs0pE5ArgCoChQw6s1mk9LxPiXpqNkn3H+SCGN7TfaRDwbICDzUap\nhIZQpO03vG4aJCILsr6eraqzu/D+l4GRqrpbRN4H/BUY193BVCK4NwDZSTM8/Vg76W9yNsDkCVO0\nAucNFK8FeLWnDVZ7i1jwfvWdYVW4J7yjqjOKPNdpBqrqzqzP54jIr0RkUCnvLaQSPe75wDgRGS0i\nUeB84MEKHLcmea0P7saFS6he3xu8OfOkkOyLmWC9cA/pNANF5AARkfTnM3Gyd2sp7y2k7OBW1QRw\nNfAosBy4O927MWXwUoBX+8KlG+EN3px5UozNSPGOYhkoIleKyJXpl30YWCIirwI3Auero1v5KarV\n71pMnjBF77r1vqqf189Wrd4XYm62UJZt2gJUZ8Xl8t3Va5vk2h5b6Zzbo+2TQlLxRe2+DkorpXfj\ncS910KYoycRDp+of7/97t98/c9zIssdQSbZy0ie8UoFXs3XixnzvDC9PHSzGWim1w4LbZ7wQ4LXU\nOvFT+yTT2E/8AAAgAElEQVSbtVKCzYLbp9wO8GrvNuhWeIM/q+9sFuLBY8Htc7kBXm1uhbdV391T\nqJViIe4/FtwBkbuhVTW50fcG96pvP0wd7EyxfriFuD9YcAeMW+2Tau806GZ4g//bJ9ksxP3HgjuA\n3Ox/10rfG4LTPslmIe4PFtwB5lb/u1amDGYEpX2Sy0Lcuyy4a4Ab1XctLJXPFaT2SS4LcW+x4K4R\nblTftdb3hmC2T3JZiLvPgrvGBL369kJ4Q3DbJ7myQ9zmiVePBXcNcqv6htrpe2cEuX1SSLEQtyCv\nLAvuGlbt6rsW+96Q3z6phQAHa6n0JAvuGlft6rtWwxtqo/9djIV4ZVlwGwAL7yqqlf53MR31xS3I\nS1O1e04a72u7fVoVbp1WC7dG64zX73tZLcXurwnB2VO80qziNnms+q6eWu1/F1OopWLyWXCbgiy8\nq8sCPF92gJv2LLhNUdXccdDC21HLFzBN6Sy4TaeqVX27Ed5eme+dq9YvYJqOWXC7LJFIsHdXM27c\ntLkrghre4N3qG2pvAY8pjQW3SxKJBPff+He+ds71fPtDP+a75/+Mxc8sd3tYHbLwdof1v00uC26X\n3POTvzHv4VdItCZIJpLs3LKTO79/L28sftPtoXWoVsLbAtx4mQW3C/buauaVJ5cQb020ezzemuDx\nO+YWfV/znmbmPfIKc+95nrffqF6g5Qp6eHu5+gYLcGMLcFyx452dhCMhkvFk3nNb179b8D2rF63l\n1mvvRFVJJVPU3SpMP20qH/7iOYhITw85T/ZinZ5eqLNs0xZWvb6pxxfpZJvYNJTlu99mxYbNnlis\nU0gmvDOLeGp1AU8tciW4W/fGWLXwjW69d+y00RUeTfUNHDqAVDL/YqQIjJgwLO/xZCLJ77/+F2LN\nsX2PAS//YzETjzmEKcdO6MnhdmjsmME9vtIye5WlhXe+fdX3SqB2V2DWEleCu74xwpjJw7v13lID\n38sBX98Q5eQLjuOfdz1LvDXe9ni4IcppF52Y9/o3lrxFKplfncdb4syb87KrwQ3p8F69mdXVqL5d\nCG/w1lL5YizAa4fvWiWlBP7qpes7DHgvhPrpF51I//378uSdT7N7+16GjRvK0DFDePR3T3LgmAM4\n+pwj6NO/CYBUqvhUwVQiVa0hdyjI4Q3+qb7B9kCpBb4L7lJ0Fu7FQr2agS4izDxzOjPPnM7Wt9/l\nhqt+w5vL15NoTbD8hdf4193PcfVNl3PAQYMZPWUEhaI70hBhxnsPq9qYO2Ph7R25/W+wAA+SmpxV\nMmby8IIfqxa+kfdRDX+96RH27momkZ5lkogladnbyr0/fQiASDTChdd9iEh9mFAk5DzWEOGQIw7m\n0BMnVWWMparmjJNqzjbJ8PqMk1w2AyWYyqq4ReR/gXOAGLAauFRVt1diYG7IrdSLtVwqXZm/tmA1\neSW1Or3tZCpJqC7EpGPGc82fPsfL/1jM3p3NjJ85ljFTR7kyo6QztVB5gz/63hlWgQdLua2Sx4Fr\nVTUhIj8ErgX+u/xheUOhlktumFcixCPRcMGpgaFwiDrZ90tRv0H9OPljx5d9vngsznMPzOflJxYR\nioQ55pwZHHHG1HbnKlfQwxv81TrJsAAPhrKCW1Ufy/ryBeDD5Q3HG1SVtcvWs+z514g2RJh+6qHs\nN3QAkB/muRV5d4L8yFmH89zf5pOM7QvvUCTE9FOn5FXUrS0x5s15mSXPrKBpYBPHf2AmoyePLPlc\nyWSSX/2/37Nxzaa2BUBvr9nE6y+v4ePXfrDLY++Ihbd3WYBXlojMAm4AQsBvVfX6nOf/A6eoFWAX\ncJWqvpp+bm36sSSQUNUZnZ2vkhcnLwP+UuxJEbkCuAJg6JChFTxtZakqf/7hX1k0dxnx1jh1dXX8\n4/a5fOTL53LEafl7A2cHeXer8TM/eQr/XrORt5atR0J1aEoZevAQzrv6fe1e19oS44YrZ7Nt03Yn\ndAWWPbuCs688g+POm1nSuZY9t5KNa7e0W7UZb4mzeO4yNl5wPAeMqmz4WHh7mwV4+UQkBPwSOB1Y\nD8wXkQdVdVnWy94ATlTVbSJyJjAbOCrr+ZNV9Z1Sz9lpcIvIE0Ch/5PXqeoD6ddcBySAO4odR1Vn\npwfL5PGTPbsV3soFq53QbnHmV6eSKVLJFP/34weZdMwhNPZuKPre7oZ4tD7KVT/5BBtWb2TT2s3s\nP2I/RhySvxDnhYde2hfaAOosk3/o5seYccZh1DfWd/r9vfbyG8SzFvJkW/Pq2ooHN+wL757mdniD\nv/re2SzAyzITWKWqawBE5M/AeUBbcKvqc1mvfwHo3kKWtE6DW1VP6+h5EfkEcDZwqnp9b9ISLHxy\nSVtoZwuFQry2YDWHnTi5pOPkhvhL/1jEqlfW0n9QE1OOm8iwcfm/dQwbcwDDxhT/x7Lk6eV5+5sA\n1IVCvLViA+MOP7jTcfUb1IdQgeX2daE6mgY2dfr+7sqssOzJqhvcDW/wd/UNwV1G39qSKHcW0iAR\nWZD19ex0MQowDMieRrWe9tV0rsuBh7O+VuAJEUkCt2Qdt6hyZ5XMAr6C8yvA3nKO5RUSEqcLlfsj\nSJyA7I7Nb2zkgZseISV1qCpP3vUsk48dz8Xf/GiXjtO7X6+Cj6dSKXo1NZZ0jCPfexj/uP1pkrQP\n7lAkxMSjD+nSeLrKwts/grYKszESbts6oZveKaX33BkRORknuLNnGRyvqhtEZDDwuIisUNXiu81R\n/jzum4A+6ZMtFJGbyzye6448YxqRaCTv8VQqxSEzOq9oc+18dxcP3PQIiViCVGsMjcVJ7G1h6XMr\nee7BeV2aK378B48iUp/zs1ag3359OHBsaf+w+g3qx2U/uIDe/XsRbYwSaYgwcOgAPvOzS4lEqrMe\nq1r3sXRjnneG17eILZXdiackG4DsamR4+rF2RGQq8FvgPFXdmnlcVTek/7sZuB+n9dKhcmeVjC3n\n/V508NSDOP6DR/H0vS8AIHUCCpd866PUN0S7fLzlL7zuHCNHormVLW9uYfDI/dvCu7OLmWOnjWbW\nZafw8K3/IBwJk0wpfQf25lM/urBL87nHHX4w37jnS2xcs5lQJMSQkftXbT54tS5Wgjcqb8D31TcE\nrwKvsPnAOBEZjRPY5wMfz36BiIwE7gMuUtXXsh7vDdSp6q7052cA3+nshIFc8l6usz51GjPPPJyV\n81YRbYwy5fgJ9OpTWisiV6HQBmfJe12dtPXCsy9mdhTgJ37kWGa+73DeXLaBXn0bGXHIgUVDd9vm\nHSx4fCG7t+1l/BEHM/6osYTqnHZPqC7EsLHuzO6ppfCGYLROMizA86XXsVwNPIozHfA2VV0qIlem\nn78Z+AawH/Cr9L/XzLS/IcD96cfCwJ2q+khn5xQ3ridOHj9Z77zl7qqf1w27d+zlex/7KYlY+4uK\nkWiYq2+6PC88Vy9dD5S/sGfFvNf5wzf/QiqlJONJIo1RRowbyhU/vohw2Bs/rzMzTXo6vJdt2gLg\naniDU3mD/2acdGR7rOcD/KCmaS+V21+ePGmq3nn737r9/mlHjCp7DJVUk3uVVFNTv1587L/fTzga\nTu81EiYcDXPKhScUrHhz903pjmQiye3fv5d4a6Jt9ki8Oca6lRuY//DCsr6fSgrynXQK8fJd5bvL\neuDu8EbpFXCHnzyFsdNGsXjuchKJJJOPHd+2ErOYMZOHl9w+ybXutX+jyfztXuOtCV56/FWOOccz\nhUN153i7cCedQoLUOsmwrWSryyruKukzoIljzzuSEz50dKehnZGpvsFZWr97xx6efWAeT9wxl7dW\n5F20bhMK11GsAxaOeu9n9dgxg6s20wTcr7zBf7sMlsJ2Iqwe7/0r9pDXXl7Dv+5+np3v7GTCzLGc\n8NFj6dO/d9XHMWbycJ66bx6/+e/bAWc15z9uf5op75nABdd+IG9zqGHjhtLQu77drc7A2Qr26LOm\nV23cXVW1i5VWefcoW4XZ86ziLuLZB+bzu+vuYuW813l7zSbm3vMCP7n8V+zavqfqY4nHEzzxu3+Q\n2NtCMp5EU0q8Nc6SZ1aw9JkVea+vkzou//7HaejTQLRXlEi901+fdtJkDjtpStXHX4pq9bvBe5V3\n0PreGVaB9xwL7gJirXEeuuWxdveDTCaS7N3VwlN/ebbq43lj8VttrQ+N7xtTvCXO/EcLX2wcNm4o\n37j7i5zwoaNpbGogHkuw/MXXmXvvC3h1Z4JaDW8IZuskwwK88iy4C9i4dgt1dfl/NKlEkpXzVrkw\novayw7ujDP73qo386+7n2Ll1Nyjs3raHR259kkf/8M8qjLJ7MuFdDRbe1WUBXjkW3AU09e9FMpF/\nYwOAvj24EVMxow8dSe4SG43HCYdgxPhhRacNPvK7f+ZtShVvjfOvu58nHsvfSMsrqnWxEiy83WAB\nXj4L7gIGDunPiPHDqAu1/+OJ1Ec48WPH9dh5VZWVC1Zzz88f4qFbHmPjWucfbyQS5qJvfZRIfYRI\nfRipqyNSH2HK8RM59fxjgMI3QN60dkvRc+16d3fPfBMVVMvhHdS+dzYL8O6zWSVFXPKdj/H7b/yF\n9Ss3EAqHSKVSnPWp0xg/Y0zZx962eTub3nyHQcMGMujAgQCkVPnTt/+PFfNeJ94SR+rqePb+eZx7\n9SyOOXsG42eM4bq7/pOFTy2leXcL42eMYeQEZ8/u7Dnf2fO9B48cVDigFfq48JtDV1RzWTx4a7ZJ\nRlBnneSyWShdZ8FdRFO/Xlx9w6W8u3E7u7fv5oDRQ4jW5+8a2BXJRJI7r7+fJc+sIBIJkYgnGXv4\naC7+1kdZ9cobrJi/qm0vcE2liMdSPHDTI0w9cRK9+/SiqX9vjn9/4Y3DCoX3rEtPZvZX/tSuXRKp\nj3DCh48uuAOi11RrcU5GJry9pFbCGyzAu8JaJZ0YeEB/Rk4YXnZoA9x3w99Z/K9lJGMJWva0kogl\nWPXKGzz460dZ+M8lBe9MUxeq4/UFa0o6fvZiHYDRhx7EJ757PoNHDgKgV99GTr/kRN572cllfy/V\nUs1+N7i/HWwhQdkitlS5LRSTzyruKnnglw/z4t9fzns8EUuw4JGFHHHGYSCSN01ERLq02jG38h4/\nYyxf+f3VqGrVtm7tCdVqmYA3dhTMFaQtYkuVCW+TzyruKnhj6TpeeCg/tDMS8SQzzjiMSLTwHXYO\nyemrx1rjvPqvpcx7+BW2b9mZ9/rcyhvwdWhXc353Nq9V3lAbs05M5yy4q+DVp5YQj+XfKzLjwLEH\nMHrKSE698ARnF8GGKNFG5+PS713Qrk3zxuI3+faHfszdP3qQ+3/xMP9z0Y089sen8o5ZKLz9rJrz\nu8F7M02yWXgba5VUQaHFPBmhcJgP/7+zADjtP07gyPdOY+X81UQbo0w8ely7u+7E4wlu/eqdtO5t\nbXeMf971LOOmH8zoKSPbPZ5pmwRJ1VsmHptpkpF90RKCtce36ZxV3FVw+KmHEonkt0Gkro7P//Jy\nhh1yIEufW8n9N/6deXNeYdwRBzPtpMl5t0pb9fIbBVdKxmMJ5s15peC5M3t7B4EbLROvV95Wfdcm\nC+4qGHHIgZx0wfGEo86NFCL1zs0U/uNrH2LI6P255Yt/4Pbv38uzf53PE7fP5UeX/ILlL76ed5yi\n7RZVYq35M1KydSW8k6kkby5fz9pl60gmC68gdUu1Wybg7fAGa53UImuVVNCW9Vt5a8UG+g/ux8GH\njmx3QfC9l5zE9NMOZfnzrxOOhjn0hIn06d+b5/+2gHUr/922oVUykSSZgDu+fy/fvu+/CIX3Verj\nDh9dcCl+pCHCtA52/Su2QKeQNxa/ye+/+RfisQQChMIhLv72Rxl7WHm3UquksWMGs6qKLRPwdtsE\namu+t7GKuyKSqRS3f+9efvLJX3Pvzx7i1mvv4IeX3MSOrbvavW7/YftxwoeP5thzZ7Tt6/3SE4va\n7UKYoUll3cp/t3ussamBD3zhfUSi4bbl+JHGKGMPP5jJx0/ocIylXKxs3t3Mb665gz3b9xLbG6N1\nb4y9O5u59dq72LNzb+d/EFVW7Vkmfqi8a2GpvLHgrojnH5jP0mdXkIgliDXHiDXH2Prvbdz+nXs6\nfW+xOdqKEirQFz/qzOl84ZYrOOHDR3P02UdwyTc/wqXfO5+6Eqb7ZcK7mFefWlZ4y1dNsfCfSzo9\nfjW5NUUwE95eZq2T4LPgroBnH5iXVzVrKsVbK9azu5MbLxx91hFEGvJXZTb0rmfYuPybCQMccNBg\nzv70GXz4i+cwYea4kkJbVXlr5b9JtsRY9uJrBV+zd+dekvH8Pno8lmDPDu9V3G70u8GbqytzWXgH\nm/W4KyDeUviioYh0un3qYSdO4rWXVvPy44sACIXqqAs5d7ApJZBz7di6izeXrqOpf29GHTqSOhG2\nrN/Kb/77T+zavpc6EZIqHPf+Izn3M7PavXfM4aMJ3R4mlWw/5kh9hHGHe6fHnauaUwSzebXfnWF9\n7+Cy4K6AQ0+YyLN/nU8q58Jhn4FN9N+/X9vX8ViCl59YxNJnV9I0sDfHnnskw8cN5aNfOpcTP3Is\nqxeupXe/RiYdM55INMzmde+w5OnlUFfH1BMmtu0kWIiq8vfZj/P0ffMIR0KoKr379eLTP7mY2V/5\nE9s27Wi3nP7ZB+Yz9aTJjJq0L/BGThjGhJljWTFv32ZXkYYI46aPZlTOHHGvqPYughlev1iZ0bZU\nfsPbgM33Dgpx4zZWk8dP1jtvubvq5+0pe3bu5edXzmb39j3EW+KEwiHqwnV88voLGTP1IMBZpn7T\n525ly/qtzratIoSjIT7whbOYOevwvGM+ccdcnvjTXFLJFOBsNnXWFafzng8eVXAMi55ezl3/c19b\n4AIgwoAh/dizY2/ejYMR4ZCjDuHT/3txu4eTqSQLn1zKvDnOEv0jZ03j8NMOJVRXeDm+V6xavdmV\nqjuzm6CXwztj+W5/hvdRg0e/pKozyjnG5ElT9c7b/9bt9087YlTZY6gkq7groHffXnz5ts/w0mOL\nWP3qGww6cCBHn3MEAwb3b3vNvDkvs2XdO21brKoq8dYE998wh8NyFttsfHMLT/xpLomsedupZIq/\nz36cyceNZ+CQfcfNePb+F9uHtnMSdmzdRThcIHRVad7dkjdFMFQX4ojTpnLEaVO7+8fhGjdaJn6p\nvMFaJ0FiFycrpL4hyrHnzuCir3+EMy8/tV1oAyz617K824iBU0mvW76h3WOLn17eVmlnS8QSPPvX\neQXP/86Gdws+HgrVkYwXmPtdH+HI0w4t+v34jVuzTMD70wSz2UXLYKhIcIvIl0RERWRQJY4XRI19\nGgs+nlKlvnd9u8cK7O7a5vkHF5BMtQ/1PTv3Fr8VmcJ7LzuZSP2+X64i9WH2O3AAM957GGAbUVWC\n38Lb5nv7W9nBLSIjgDOAt8ofTnAde96R+dP+xLnTzvCcaX+HnTi56HEUZc2ra9s99vpLa4rOBx85\n4UBOOf94Lv+fCzn0hEmMmTaKsz59Op//1aeI1kc6ndvtR25U3eCPOd7ZrPr2r0pU3D8DvgJU/yqn\nj4yfMYZTLjieUCRMfa96oo1R+g3qy6d+eGHeXtn7D9+PwSP3K3icOqmjZXf73QFD4RCSdx94QIT9\nD3LCZOy0UVzyrY9y1U8/wfHvP6rdVrG2EVXl+GGOdzYLb38q6+KkiJwHbFDVVzvbqF9ErgCuABg6\npPDCkqA7/aITOeacGbyx+C169evF6PQ860JOueA9/N9PHmx3gRKcvUwOnnZQu8fGzxyLFvi5GYmG\nmXHGtJLHV8peJn5Q7XtVFuKHi5UZtkWs/3RacYvIEyKypMDHecBXgW+UciJVna2qM1R1xoB+A8od\nt2819e/Noe+ZyJipB3W4wGbayZMZPm7ovvaKCJH6MGd+8lR69+nV7rXR+giXfPtjRBoiRBujROoj\nhKNhTrngOEZNKq0VYi2TyvFTvzvDtoj1l04rblU9rdDjInIoMBrIVNvDgZdFZKaq2h0+yxQKh7jy\nZ5ew8MmlLJq7jMamBo45ZwajJhee7jZ+xhi+ec+XWfrcSmItMSbMHJs3s6UUVnVXhp+mCWazKYPd\nIyKzgBuAEPBbVb0+5/kJwO+A6cB1qvrjUt9bSLdbJaq6GGj7Pysia4EZqvpOd48ZBMlUipXzVrFu\nxb8ZMKQfh500ifrG+s7fWEA4HGbGGYcx44zDSnp9Q6/6suZfB+2OOW5s/5rNwrs2iEgI+CVwOrAe\nmC8iD6rqsqyXvQt8Hnh/N96bx+ZxV1BrS4wbP/Mbbv/uPTz+x6e4/xdz+N7HfsbGtf761TMoFyoz\n3GqZgD/bJmBTBrtoJrBKVdeoagz4M3Be9gtUdbOqzgdyNy/q9L2FVGzlpKqOqtSx/OrJO55m09ot\nbRcU4y1x4hLnju/dy5d+e5XLoytNIKtuly9UZipvPwpK9d3aHCv3B/ggEVmQ9fVsVZ2d/nwYkH3w\n9UDhvSnydeu9tuS9gl5+YnHeLBAUNq/byq5tu+kzoMmdgXVDUHrdGW7tIJgxacj+LPNZyyQjCOFd\nXx8ud4HWO17aq8RaJZVUdJKIOsshfSJoM0zcntudzW8tkwybcdKhDUB2VTA8/ViPvdeCu4JmvHda\n/gpGEYYcNLjtVmXGHW4uh8/wa787w/reRc0HxonIaBGJAucDD/bkey24K+jk849j2NgDiDZGkbo6\noo1Revdt5MKvf8jtoXVZkFZTZnO76vZ7eINV37lUNQFcDTwKLAfuVtWlInKliFwJICIHiMh64IvA\n10RkvYj0Lfbezs5pPe4KitZHuPoXl/P6K2+wfsUG+g/px5TjJ7ZbXm7c44ULleDfaYLZgtD3riRV\nnQPMyXns5qzPN+K0QUp6b2es4q4wEeGQ6Qdzysffw/RTp/o+tINWdY8dM9j1qhus8jblseA2RQXt\nImU2L4W3n1nf2x0W3KZTQay6vcJvuwkWY9V3dVlwmw5Z1V0dFt6mKyy4TU3yWtUNFt6mdBbcplNB\nnRrolQuVELzwtr53z7LgNjXPwrtnWPXdcyy4TcmCWnV7SRBmmmSz8O4ZFtymJEG+SAneqbohODNN\nMrLD2wK8Miy4Tc3zWtWdEbTwtuq7ciy4TZcEsV0C3rpQCcHrd2dYeFeGBbcpWdDbJeC9lglYeJt8\nFtzGpHmxZRLk8LYpg91nwW26LKjtkgwvVd0QvJkm2az67h4LbtMlQW+XeLHqhuDNNMlm4d11FtzG\nFOC1qjvDwtuABbfppiC3S7xcdUOww9v63qWx4DZdFvR2CXhvemBG0MMbrPouhQW3MR2w8HaHhXfH\nLLhNtwW5XQLebZlAsGeaZGSvtjTtWXCbbqmFdkmGF6tuCPZME9MxC25jOuDlqjvDwrv2WHAbUwIv\nV91g4V1rLLhNtwX1zji5vF51W3jXnrKDW0Q+JyIrRGSpiPyoEoMyxmu8Oj0ww8K7tpQV3CJyMnAe\ncJiqTgZ+XJFRGeNRfghvE3zlVtxXAderaiuAqtqkyxpUC+0S8H7LBGymSa0oN7gPAd4jIi+KyL9E\n5MhKDMr4Ry1NC8zwctWdYeEdbJ0Gt4g8ISJLCnycB4SBgcDRwH8Bd4uIFDnOFSKyQEQWbNuxraLf\nhDHV4peqGyy8g6zT4FbV01R1SoGPB4D1wH3qmAekgEFFjjNbVWeo6owB/QZU9rswpsq8XnVbeAdb\nua2SvwInA4jIIUAUeKfcQRl/qZVpgRl+qLrBwjvIyg3u24CDRWQJ8GfgElXV8odljLd5fXpghs00\nCaaygltVY6p6Ybp1Ml1Vn6zUwIwxlWEzTXqeiMwSkZUiskpErinwvIjIjennF4nI9KznbhORzekC\nuCS2ctKYbvJL1Z1h4d0zRCQE/BI4E5gEXCAik3JediYwLv1xBfDrrOd+D8zqyjktuE1F1FqfO5sf\nwtv63T1qJrBKVdeoagynbXxezmvOA/6YnsjxAtBfRIYCqOpc4N2unDBcgUEbU7PGjhnMqtX+WHc2\nacj+LNu0hVWvb2LsuCFuD6eqWvfGyi0sBonIgqyvZ6vq7PTnw4Dsn97rgaNy3l/oNcOAt7szGAtu\nYypg9ZJ1jJkywu1hdCoT3rWmvjFS7mKxd1R1RqXGUy5rlRhTJr9MD8ywi5UVtwHI/qk9PP1YV19T\nMgtuUzG13OcGf/S6s1l4V8x8YJyIjBaRKHA+8GDOax4ELk7PLjka2KGq3WqTgAW3MRXhx6obLLwr\nQVUTwNXAo8By4G5VXSoiV4rIlemXzQHWAKuA3wCfybxfRO4CngfGi8h6Ebm8s3Naj9uYCvJLrxtq\n+2JlpanqHJxwzn7s5qzPFfhskfde0NXzWcVtTIX4reoGq7z9yoLbVFwt97n9tigHbFm8H1lwm4qq\nxf25C/FjeFvV7R8W3MZUmB9bJhkW3v5gwW1MD/Fj1Q0W3n5gwW1MD/Br1W39bn9wZTrgsteWvTPt\n5ClvVuhwg/DmzRu8Oi7w7thsXF3j1XFBZcd2UIWOExiuBLeqVuzHuogs8NIeAhleHRd4d2w2rq7x\n6rjA22MLAmuVGGOMz1hwG2OMzwQhuGd3/hJXeHVc4N2x2bi6xqvjAm+PzffE7u1rjAm6yeMn6523\n3N3t9087ecpLXurZB6HiNsaYmhKY4BaRz4nIChFZKiI/cns82UTkSyKiIjLI7bEAiMj/pv+sFonI\n/SLS3+XxdHiHbLeIyAgR+aeILEv/vfqC22PKJiIhEXlFRB5yeywZItJfRO5J//1aLiLHuD2mIApE\ncIvIyTg34zxMVScDP3Z5SG1EZARwBvCW22PJ8jgwRVWnAq8B17o1kBLvkO2WBPAlVZ0EHA181kNj\nA/gCzv7PXnID8IiqTgAOw3vjC4RABDdwFXC9qrYCqKqX7t76M+ArgGcuJqjqY+nN3wFewLmNkltK\nuUO2K1T1bVV9Of35LpwQGubuqBwiMhw4C/it22PJEJF+wAnArQCqGlPV7e6OKpiCEtyHAO8RkRdF\n5AlZqNUAAAMrSURBVF8icqTbAwIQkfOADar6qttj6cBlwMMunr/Y3a89RURGAYcDL7o7kjY/xykI\nUm4PJMtoYAvwu3QL57ci0tvtQQWRb+6AIyJPAAcUeOo6nO9jIM6vs0cCd4vIwVqFKTOdjOurOG2S\nqutoXKr6QPo11+G0A+6o5tj8RkSagHuB/1TVnR4Yz9nAZlV9SUROcns8WcLAdOBzqvqiiNwAXAN8\n3d1hBY9vgltVTyv2nIhcBdyXDup5IpLC2Sthi1vjEpFDcSqQV0UEnHbEyyIyU1U3ujWurPF9Ajgb\nOLUaP+A6UNG7X1eaiERwQvsOVb3P7fGkHQecKyLvAxqAviJyu6pe6PK41gPrVTXzW8k9OMFtKiwo\nrZK/AicDiMghQBSXN99R1cWqOlhVR6nqKJy/1NOrEdqdEZFZOL9mn6uqe10eTil3yHaFOD9xbwWW\nq+pP3R5Phqpeq6rD03+vzgee9EBok/67vU5ExqcfOhVY5uKQAss3FXcnbgNuE5ElQAy4xOUq0utu\nAuqBx9O/Dbygqld2/JaeoaoJEcncITsE3KaqS90YSwHHARcBi0VkYfqxr6ZvDGsK+xxwR/qH8Brg\nUpfHE0iBCO70bATXK46OpKsjT1DVsW6PIVuhO2R7gao+A4jb4+iIqj4FPOXyMNqo6kLAMysMgyoo\nrRJjjKkZFtzGGOMzFtzGGOMzFtzGGOMzFtzGGOMzFtzGGOMzFtzGGOMzFtzGGOMzFtzGGOMzFtzG\nGFOmzu7iJI4b088vEpHppb63EAtuY4wpQ4l3cToTGJf+uAL4dRfem8eC2xhjylPKXZzOA/6ojheA\n/iIytMT35rHgNsaY8pRyF6dir+nWHaACsTugMcZ0ZNlryx6ddvKUQWUcokFEFmR9PVtVZ5c7ru6y\n4DbGBJ6qzurBw5dyF6dir4mU8N481ioxxpjylHIXpweBi9OzS44Gdqjq2yW+N49V3MYYU4Zid3ES\nkSvTz9+Mc6OQ9wGrgL2k7wzU3TtAid3hyxhj/MVaJcYY4zMW3MYY4zMW3MYY4zMW3MYY4zMW3MYY\n4zMW3MYY4zMW3MYY4zMW3MYY4zP/H9p0ibOMD6UjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff478c0fe80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "vlr = VariationalLogisticRegressor()\n",
    "vlr.fit(X_train, y_train)\n",
    "y = vlr.proba(X).reshape(100, 100)\n",
    "plt.figure(figsize=(6,6))\n",
    "plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)\n",
    "plt.contourf(x0, x1, y, np.array([0., 0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99, 1.]), alpha=0.25)\n",
    "plt.colorbar()\n",
    "plt.xlim(-7, 7)\n",
    "plt.ylim(-7, 7)\n",
    "plt.gca().set_aspect('equal', adjustable='box')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
